%\documentclass[5p, authoryear]{elsarticle}

% A utiliser pour la soumission :
\documentclass[preprint, review, authoryear, 12pt]{elsarticle}

\usepackage{rotating}

\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}
\usepackage{lmodern}

\usepackage{amsmath,amsfonts,amssymb}
\usepackage{textcomp} % Symbole celcius, euro etc.
\usepackage{subfig}
\usepackage[babel=true]{csquotes}
\usepackage{array}
\usepackage{multirow}
\usepackage{url}
%\usepackage{hyperref}

\graphicspath{{img/}}

\newcolumntype{C}[1]{>{\centering}m{#1}}

\journal{Computers and Electronics in Agriculture}

\usepackage{lineno}
\linenumbers

\selectlanguage{english}

\begin{document}

\begin{frontmatter}

\author[INRA]{Jean-Romain Roussel}
\ead{jromain.roussel@gmail.com}

\author[INRA,AgroParisTech]{Frédéric Mothe}
\ead{mothe@nancy.inra.fr}

\author[Loria]{Adrien Krähenbühl}
\ead{adrien.krahenbuhl@loria.fr}

\author[Loria]{Bertrand Kerautret}
\ead{bertrand.kerautret@loria.fr}

\author[Loria]{Isabelle Debled-Rennesson}
\ead{Isabelle.Debled-Rennesson@loria.fr}

\author[INRA,AgroParisTech]{Fleur Longuetaud\corref{cor}}
\ead{longueta@nancy.inra.fr}

\cortext[cor]{Corresponding author}

\address[INRA]{INRA, UMR1092 LERFoB, 54280 Champenoux, France.}
\address[AgroParisTech]{AgroParisTech, UMR1092 LERFoB, 54000 Nancy, France.}
\address[Loria]{LORIA, UMR CNRS 7503, Université de Lorraine, 54506 Vandœuvre-lès-Nancy, France.}

\title{Automatic knot segmentation in CT images of wet softwood logs using a tangential approach}

\frontmatter

\begin{abstract}
Computed Tomography (CT) is more and more used in forestry science and wood industry to explore internal tree stem structure in a non-destructive way. Automatic knot detection and segmentation in the presence of wet areas like sapwood for softwood species is a recurrent problem in the literature. This article describes an algorithm named \textit{TEKA} able to segment knots even into sapwood and other wet areas by using parallel tangential slices into the log that enable to follow the knot from the stem pith to the bark.
On each tangential slice, knot pith is detected, then knot diameter is estimated by analyzing gray level variations around the knot pith.
A validation was performed on 125 knots from five softwood species. The CT slice resolution ranged from 0.4 to 0.8 mm/pixel with an interval between slices of 1.25 mm. Compared to manual diameter measurements performed on the same CT slices, the \textit{TEKA} algorithm led to a RMSE of 3.37 mm and a bias of 0.81 mm, which is rather good compared to other algorithms working only in heartwood.
\end{abstract}

\begin{keyword}
Computed tomography \sep Sapwood \sep Knottiness \sep Algorithm \sep Wood quality
\end{keyword}

\end{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
\section{Introduction}

The application of X-ray Computed Tomography (CT) to wood science and industry was investigated since many years. After the pioneer works of \cite{Taylor1984}, \cite{Funt1985} and \cite{Funt1987}, numerous works have been conducted to develop algorithms to automatically detect wood features in CT images with a special interest for knots. \citet{Fleur01} presented a review of the literature about automatic knot detection algorithms. Recent works can be added to this review: \citet{Aguilera2012} (species not specified), \citet{Breining2012} on Norway spruce and \citet{johansson2013} for Scots pine and Norway spruce.
In overall, algorithms for knot detection and segmentation are efficient on dried wood, but a recurrent problem mentioned in the literature is the presence of wet areas generally within sapwood. 
For many softwood species, sapwood has a much higher wood density than heartwood at fresh state \citep{Polge1964} due to the higher water content of sapwood. Within fresh logs of these species, classical approaches based on gray level values are not efficient because fresh sapwood has almost the same density as knots.\\

\citet{Aguilera2012} is the continuation of an approach based on simulated annealing in deformable contours \citep{Aguilera2008b, Aguilera2008a}. Using deformable contours for knot segmentation is an original approach that can work in presence of sapwood. However, in the examples of their experiments, only a very small part of the knots is included within sapwood. Moreover, the segmentation process is not fully automatic since the deformable model must be manually initialized and the method was not statistically validated.\\

\citet{Breining2012} algorithm is a classical approach based on gray level thresholding. They first remove sapwood in the CT images in order to detect knots in heartwood only, based on a fixed gray level threshold corresponding to a density of 900 kg.m$^{-3}$. Morphological operations are then used to improve the knot detection. The algorithm was designed to work within fresh heartwood but not within sapwood. Indeed, the method shows some weaknesses in presence of partly dried sapwood. A statistical validation was performed based on 55 knots from 55 cross-sections. The knots with very unclear border were avoided in the validation sample, which could artificially improve the accuracy results.\\

Until now, \citet{johansson2013} are the only ones to propose a method designed for working in heartwood and in sapwood as well. Their algorithm is the continuation of \citet{Grundberg01} works. The algorithm is also based on images of concentric surfaces (CS) or cylindrical shells within the logs, following approximately annual rings. The main difference is that \citet{johansson2013} have applied their algorithm on images of lower resolution than classical CT images obtained from medical scanners. The objective was to process images like the ones which would be obtained by a high speed industrial CT scanner. The knot detection is based on 10 CS with a minimum of five CS in heartwood. CS are thresholded in order to detect high density objects and then ellipses are fitted on these objects. Ellipses which can be matched through at minimum three consecutive heartwood CS are assumed to correspond to knots. Regression models for size and location of knots are fitted from the detections in heartwood and they are then used to generate sub-images in sapwood CS supposed to contain the knots. Computation of gray level standard deviations in rows and columns in these sub-images confirm or not the presence of a knot. If a knot is present they \enquote{try to find the position and size of it in the sub image using morphological dilation}. This last step is not detailed. Since the authors write that the detection in sapwood succeeds only for \enquote{knots that have higher density than the surrounding sapwood}, it may be supposed that knots are detected and measured by gray level thresholding. Given our images, a method based on a threshold could fail in many cases.\\

This paper presents an algorithm - named \textit{TEKA} - designed for knot segmentation into wet logs in which sapwood can have a density similar to knot density. We chose to focus on the segmentation step.
That means that we started from already defined \enquote{knot areas}, i.e., angular sectors, radius range and slice interval framing each knot.
\textit{TEKA} is then able to separate the knot from sapwood or moisture areas within each knot area. Indeed, in our approach, we have assumed that detection (i.e., localization of knot areas) and segmentation (i.e., segmentation of the knot within each knot area) were two different steps.
The algorithm could plug after any existing software able to achieve the detection step and to define knot areas.

Starting from a knot area, \textit{TEKA} uses an original approach by looking at the log in a tangential view rather than the classical transversal view (i.e., CT slices or cross-sections). A tangential image is orthogonal to the CT slices and tangential with regard to annuals rings. \citet{Grundberg01} followed by \citet{johansson2013} already presented a quite similar approach but based on concentric surfaces centered on the log pith (see above).
The main difference is that our segmentation method is based on knot pith detection and analysis of gray level variations around the knot pith rather than on classical image thresholding.
Furthermore, it was designed for working as well in heartwood and in sapwood. Validation results obtained on five softwood species are provided and discussed hereafter.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
\section{Materials and methods}

\subsection{Sampling}

We developed and validated the algorithm based on 125 knots from 16 wet logs (12 trees) of five softwood species: four Douglas-fir logs (25 knots), three silver fir logs (29 knots), three European larch logs (26 knots), three Scots pine logs (20 knots) and three Norway spruce logs (25 knots). Logs were provided by Siat Braun sawmill (France).

The diameters of the logs ranged from 12 to 27 cm with a sapwood ratio ranging from 24 to 62 \% of the log radius.

Before being transported to the laboratory, the logs were stored during 1 to 5 months under water sprinkling on the log yard. In consequence, the moisture content of logs could be variable: from fresh, with an entire sapwood visible in CT images, to partly air-dried, with a sapwood that could have partly dried in some areas (Fig. \ref{imagesBillons}). The algorithm was designed to be able to process all the range of drying states which may be encountered in a sawmill.

\begin{center}
*****Figure \ref{imagesBillons} about here*****
\end{center}

The knots have been selected manually trying to be representative of the usual sizes for each species. The range of knot diameters (maximal diameter of the knot along its radial profile) of our sample was 5.6 mm to 42 mm with an average of 19 mm. We selected only knots originating from the pith and ending close to the bark.
Figure \ref{Images_zone_noire} shows transverse, radial and tangential cross-sections of some of the sampled knots.

\begin{center}
*****Figure \ref{Images_zone_noire} about here*****
\end{center}

\subsection{CT scanning}
The logs were analyzed using a medical CT scanner (BrightSpeed Excel by GE Healthcare). Stacks of 512 $\times$ 512 pixels images were obtained. Image thickness and interval between two images were 1.25 mm. The pixel width ($w$) ranged between 0.4 and 0.8 mm/pixel depending on the log diameter.

\subsection{Manual detection and measurement of knots}

For the validation of our algorithm, each knot was manually described with the \textit{Gourmands} plug-in \citep{Colin01} for ImageJ software \citep{Schneider2012}. This tool allows the user to browse the stack of images and place markers on the knot borders. One line of markers have to be placed on both sides of each knot. The knot diameter in the direction tangential to the annual growth rings of the log is assessed by computing the distance between the two lines. The middle line defines the knot trajectory. For each knot, 11 values of knot local diameters and trajectory points were recorded every 10\% of the knot length from the log pith to the bark. The maximum of the 11 local diameters of each knot was also recorded.

Here-after, it will be assumed that the trajectory points correspond to the knot pith location. This is probably correct in the tangential direction but less true in the vertical direction (i.e., along the main log axis) if the pith is not vertically centered.

A repeatability test based on 44 randomly chosen knots was performed to evaluate the accuracy of the human measurements. The knots were measured twice by the same operator.
The root mean square deviations between both set of measurements were 1.4 mm for the maximal diameter, 2.1 mm for the local diameter, 5.1 mm for the vertical pith coordinate and 0.7 mm for the horizontal pith coordinate.

\subsection{Automatic knot segmentation algorithm}

The \textit{TEKA} algorithm, written in Java language as a plug-in for ImageJ software, works with stacks of tangential images (i.e., images sliced tangentially to annual growth rings). We have resliced the original stack of CT images (cross-section slices orthogonal to the main log axis) to produce stacks of tangential images of the knots from the log pith to the bark (Fig. \ref{decoupeTan}). This step was performed using an ImageJ macro taking as input a radial line passing through the knot, manually drawn by the operator.
The tangential images were produced with a pixel size $w$ equal to the pixel size of the original CT slices. The distance between two tangential images was also fixed to $w$.

The \textit{TEKA} plug-in delivers knot radius and knot pith coordinates for each tangential image.
Segmented images of the knot are also delivered to visualize the results.

%The \textit{TEKA} algorithm is described section \ref{algo}.

\begin{center}
*****Figure \ref{decoupeTan} about here*****
\end{center}

\subsection{Statistical validation}

We have compared the automatic diameter measurements (based on the radii provided by \textit{TEKA}) and the automatic knot pith detections to the diameters and knot centers measured manually. The comparison was made for each knot at 11 positions by interpolating the data to get manual and automatic measurements every 10\% of the knot length.
In total, since the first knot diameter is null, 1250 diameters and 1375 coordinates were compared (125 knots $\times$ 10 or 11 measurements per knot).

For allowing to compare the accuracy of the algorithm in heartwood and sapwood, the stacks of tangential images were reviewed to decide for each knot at which image occurred the transition from heartwood to sapwood. In total, 589 measurements over 1250 (47\%) were attributed to sapwood.

The statistical software \emph{R} \citep{R} was used for statistical validation. The statistical values that were computed are: root mean square deviation (RMSD) and error (RMSE), mean absolute error (MAE), r-square ($R^2$) and the mean bias (computed as the mean of (automatically measured values $-$ manually measured values)).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
\section{Description of the segmentation algorithm}
 \label{algo}

The \textit{TEKA} algorithm can be divided into three main steps:

\begin{enumerate}[Step 1:]
\item Knot pith detection;
\item Knot diameter measurement;
\item Post-processing.
\end{enumerate}

Step 1 and 2 are processed on each tangential image (after reslicing the original CT images, see Fig. \ref{decoupeTan}) whereas step 3 concerns the whole profile.

\subsection{Step 1: Knot pith detection}

The \textit{PithExtract} algorithm initially presented by \citet{Fleur02} was used. This algorithm was recently improved by \citet{Boukadida01} and validated on a big set of 100451 CT-images, with pixel size ranging from 0.2 to 1 mm. \textit{PithExtract} is based on a Sobel edge detection, where edges correspond here to the border of the knot cross-section, and the Hough accumulation principle. The pith detection is robust even with partial information, noise or ellipticity.

%From the knot pith profile, the algorithm calculates an inclination profile of the knot. The inclination profile is the discrete derivative of the vertical pith position profile, assuming that horizontal deviations are negligible. In the following, $\alpha$ is the vertical inclination of a knot at a given location between the stem pith and the bark.

\subsection{Step 2: Knot diameter measurement}

\textit{TEKA} computes the knot radius on each tangential slice in two sub-steps:

\subsubsection{Sub-step 1: Polar elliptic transformation centered on the knot pith}

%Since knot cross-sections in planes oriented perpendicularly to the knot pith profile are almost circular and due to the knot inclination, it results that knot intersections with tangential images have an elliptical shape. For each tangential image, an ellipticity rate $\tau$ is computed from the $\alpha$ angle.
%\[ \tau = \frac{R_2}{R_1} = \cos \alpha \]
%Where $R_1$ is the major radius and $R_2$ is the minor radius of an ellipse.

It was assumed that knots have circular shapes on cross-sections oriented orthogonally to the knot pith profile. 
Due to the knot inclination, the knot section on a tangential plane perpendicular to the stem cross section has an elliptical shape.
For each tangential image, an ellipticity rate $\tau$ was computed from the local vertical inclination angle $\alpha$. This angle $\alpha$ is the discrete derivative of the vertical pith position, assuming that azimuthal deviations are negligible:\\
\[ \alpha = \arctan{ \frac{\Delta z}{\Delta r}}\]
\[ \tau = \frac{R_2}{R_1} = \cos \alpha \]

Where $\Delta r$ is the horizontal distance between two successive tangential slices (which was here equal to the pixel width $w$), $\Delta z$ is the corresponding vertical deviation of the knot pith, $R_1$ and $R_2$ are the major and minor radii of an ellipse.

The tangential images were first smoothed with a Gaussian blur filter with a radius of $l$ pixels\footnote{The values of the setting parameters used for testing the algorithm are given in section \ref{Choice_of_the_parameters}.} in order to smooth the gray level profile used in section \ref{substep2}.

Then, a polar elliptic transformation centered on the automatically detected knot pith was performed. Like polar circular transformation, polar elliptic transformation consists in converting the image from cartesian to polar coordinates but by using parametric equations of ellipses rather than circles. The resolution of the resulting polar image was equal to one degree in the angular direction and to $w$  in the radial direction, where $w$ is the pixel width on the tangential slices\footnote{Actually, the radial resolution of the polar elliptic images varied between $w$ and $\frac{w}{\tau}$  depending on the angular position.}.
Figure \ref{polar} illustrates how that elliptic transformation (Fig. \ref{polarc}) is more appropriate than circular transformation (Fig. \ref{polarb}) to obtain a vertical pattern corresponding to the knot (on the left of the images).

\begin{center}
*****Figure \ref{polar} about here*****
\end{center}

\subsubsection{Sub-step 2: Analysis of the gray level profile} \label{substep2}

The gray level profile corresponding to the mean values of each pixel column was computed (Fig. \ref{histogramme}). This profile $f$ has two characteristics:

\begin{enumerate}
\item A maximum value into the knot because of knot sapwood density at point $A(r_{A}, y_{A}=f(r_{A}))$;
\item A negative derivative from point $A$ to the end of the knot/wood transition that occurs at point $B(r_{B}, y_{B}=f(r_{B}))$, where the derivative is 0.
\end{enumerate}

\begin{center}
*****Figure \ref{histogramme} about here*****
\end{center}

Points $A$ and $B$ are automatically detected. $r_{A}$ is an underestimation of the ellipse minor radius whereas $r_{B}$ is an overestimation. The real minor radius is estimated by equation (\ref{eqrayon}) at point $C(r_C, y_C=f(r_C))$ (Fig. \ref{histogramme}).

\begin{equation}
r_{C} = f^{-1}\left( y_B + \beta (y_{A}-y_B) \right)
\label{eqrayon}
\end{equation}

$f^{-1}$ is the reciprocal function of $f$ reduced to $[r_{A}, r_{B}]$ and $\beta$ a parameter included in $[0,1]$.
%Using the reciprocal function allows to take account of $f$ slope. 

The radius at point $C$ corresponds to the minor radius of the ellipse. The major radius is computed using the ellipticity rate $\tau$.

Figure \ref{segmentation} illustrates two examples of knot segmentation into sapwood and into heartwood.

\begin{center}
*****Figure \ref{segmentation} about here*****
\end{center}

\subsection{Step 3: Post-processing}

The post-processing step was designed to improve the algorithm accuracy. First, outlier radii are identified and replaced by linear interpolation. Then, the profile of knot radii is smoothed. Last, the very end of the knot profile is corrected.

\subsubsection{Detection and correction of outliers}
\label{correction}

Some errors might appear during the knot radius computation when the pattern previously described was not present in the gray level profile. In such a case, a large overestimation or underestimation of the radius was observed. For this reason, the algorithm processed to a correction of the radii identified as outliers.\\

To identify outliers we used a non parametric regression method. We chose the LOWESS algorithm (LOcally WEighted Scatterplot Smoothing) created by \citet{Cleveland1979} which is considered to be resistant to outliers. We fitted the LOWESS curve on the knot radius profile (Fig. \ref{lowess}a) and computed the residuals between model and data. The outliers are the points out of the limits given by $[Q_1 - k \times IQR, Q_3 + k \times IQR]$ where $Q_1$ and $Q_3$ are the first and third quartiles respectively and $IQR$ the inter-quartile distance. The constant $k$ was fixed to 1 by trial and error method.\\

\begin{center}
*****Figure \ref{lowess} about here*****
\end{center}

Then, outliers were removed and the gaps were filled by linear interpolation between the boundaries of the gaps (Fig. \ref{lowess}b).\\

We used the standard Java implementation of the LOWESS regression algorithm from the Common math 3.2 API\footnote{\url{http://commons.apache.org}}.

\subsubsection{Smoothing of the knot radii profile}

The knot radius at the log pith was forced to be 0. Then, the profile was smoothed using an approximated Gaussian blur filter with $s$ pixels radius. We chose the approximated Gaussian filter because it is simple and it better preserves local variations than a not weighted moving average. This step allows to get a continuous profile (Fig. \ref{lowess}c).

\subsubsection{Extrapolation of the knot radius at the bark side}

The knot segmentation was often difficult in tangential slices very close to the bark. For this reason, the radii estimated on the last $p$ percent of slices located at the bark side were deleted.
The knot end was extrapolated with a constant radius value equal to the last valid value, still centered on the location given by the \emph{pithExtract} algorithm (Fig. \ref{lowess}c).

\subsection{Choice of the parameters}
\label{Choice_of_the_parameters}

For testing the accuracy of the \textit{TEKA} algorithm in the next section, the parameters were set empirically by trial and error method to the following values:

\begin{itemize}
\item The radius $l$ of the Gaussian blur mask applied to the image before polar elliptic transformation was set to 7 pixels;
\item The $\beta$ parameter was set to 0.75;
\item The percentage $p$ of the knot length for which diameters were extrapolated was set to 10\%;
\item The radius $s$ of the Gaussian blur mask for the knot radius profile was set to 22 pixels;
\item The LOWESS regression algorithm from the Common math 3.2 API takes 3 parameters named \textit{bandwidth, robustness} and \textit{accuracy} which were set respectively to 0.33, 3 and 0;
\item The $k$ constant used to define the outlier diameters was set to 1;
\item The \textit{PithExtract} parameters were set to the default values described in \citet{Boukadida01} except the wood / background threshold ($B$) which was set to -300 HU rather than -700 HU.
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
\section{Results} \label{accuracy}

Figures \ref{result1}, \ref{result2} and \ref{result3} show the segmentation results by \textit{TEKA} for three knots of various size and species. The three examples are visually satisfactory in comparison with what would be obtained by manual segmentation.

\begin{center}
*****Figure \ref{result1} about here*****
\end{center}

\begin{center}
*****Figure \ref{result2} about here*****
\end{center}

\begin{center}
*****Figure \ref{result3} about here*****
\end{center}

Table \ref{tabstattout} summarizes the results of the comparison with manual measurements by species and wood compartments for the maximal and local diameters and for knot pith positioning.

\begin{center}
*****Table \ref{tabstattout} about here*****
\end{center}

\subsection{Accuracy of the diameter measurements}

Figure \ref{diameter_plot} shows the plot of local diameters automatically measured by \textit{TEKA} versus corresponding manual measurements. All the 1375 measurements are represented. By definition the knot diameter at the log pith location was always 0 and thus 125 points are at (0,0) in the plot. Statistics were calculated by removing these points which would artificially improve the results.\\

\begin{center}
*****Figure \ref{diameter_plot} about here*****
\end{center}

For local diameter measurements, RMSE was 3.37 mm, mean absolute error was 2.26 mm, mean bias was +0.81 mm and $R^2$ was 0.85.
Figure \ref{diameter_error_boxplot} presents the errors as a function of the position along the knot. Each box is made of 125 measurements from the 125 knots. The accuracy of the local diameter measurements was almost the same everywhere in the knot although errors were slightly higher close to the bark.

% Rajouter ici quelque chose sur les outliers et la forme en U
% The majority of outliers are found close to the pith. The density of probability for outliers have a pattern like a “U” with a maximum at pith position and an other maximum à $p$\%. Over $p$\% it is the extended part and outlier notion does not make sense. On average 14\% of measures are outliers. Fig. \ref{lowess} have 12\% of outliers respectively positioned close to the pith and close to the $p$\% position so it is representative of the sample. Close to pith outliers are often the result of the few quantity of information available but in the general case the local breakage of the pattern is not easily explicable. Three main cases observed could be described (i) it could be due to an internal knot dark artefact (under estimation) (ii) it could be due to the a local absence of pattern due to a local absence of dark structure at the knot periphery (over estimation) and (iii) some outliers appear for very small variation which are not visually outliers when the raw diameter profile is perfectly followed by the lowess curve. In this case the interquartile interval is very small and very little deviation are corrected. In any case each group of outliers must be interpreted individually.

\begin{center}
*****Figure \ref{diameter_error_boxplot} about here*****
\end{center}

For maximal diameter, RMSE was 3.83 mm, mean absolute error was 2.62 mm, mean bias was +0.69 mm and $R^2$ was 0.85.

The largest errors were obtained for Scots pine (maximal diameter) and Douglas fir (local diameter), the smallest errors for larch and silver fir. No difference was observed between sapwood and heartwood.

\subsection{Accuracy of the knot pith position}

Figure \ref{moellez} shows the absolute error made on vertical and horizontal pith positioning versus the relative position along the knot. Each box is made of 125 measurements from the 125 knots.\\

\begin{center}
*****Figure \ref{moellez} about here*****
\end{center}

On vertical positioning, RMSE was 4.21 mm, mean absolute error was 2.04 mm and mean bias was -0.11 mm. $R^2$ between automatic detection and manual positioning was 0.96. For the positions 0\% and 100\%, where the errors were the highest, the RMSE were 9.44 mm and 5.38 mm, respectively. Between positions 10\% and 90\%, the RMSE ranged between 1.36 mm and 5.96 mm.

On horizontal positioning, RMSE was 1.59 mm, mean absolute error was 0.81 mm and mean bias was +0.05 mm. $R^2$ between automatic detection and manual positioning was 0.97. For the positions 0\% and 100\%, where the errors were the highest, the RMSE were 3.31 mm and 2.49 mm, respectively. Between positions 10\% and 90\%, the RMSE ranged between 0.47 mm and 1.94 mm.

The largest errors were obtained for Scots pine, the smallest errors for silver fir (vertical position) and Douglas fir (horizontal position). No difference was observed between sapwood and heartwood.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
\section{Discussion} \label{discussion}

Regarding the accuracy of our algorithm, the observed errors on local knot diameter (3.37 mm) and pith position (4.21 mm vertically and 1.59 mm horizontally) were in the same order as the root mean square differences between two repetitions of manual measurements (2.13 mm for diameter, 5.12 mm vertically, 0.73 mm horizontally). Considering moreover that the knot pith was not really measured manually but estimated through the center line of the knot, this result is really satisfactory.

The pith positioning errors were bigger vertically than horizontally. It is probably due to the uncertainty of the manual measurements which is also bigger vertically than horizontally. It can also be related to the voxel size of the initial CT images (1.25 mm vertically, 0.36 to 0.81 mm horizontally) even if no effect of voxel horizontal size was observed on RMSE.

The pith positioning errors were bigger close to the log pith and close to the bark. Close to the log pith the knots are very small, with a fuzzy shape and very few edge pixels to be used by \textit{PithExtract} algorithm to find the pith. Furthermore, other knots of the same whorl might appear in the images and lead to detection errors. At the other side, the external knot end is often perturbed by bark, resin pockets or other complications which increased the detection errors. Moreover, some of the tested knots which ended just before the bark led to incorrect results.

For comparison purposes, \citet{Breining2012} gave a RMSE of 4 mm, a bias of 1.7 mm and a $R^2$ of 0.68 based on 119 knot measurements in heartwood only and for Norway spruce. \citet{Fleur01} provided validation results based on 365 knots from Norway spruce and silver fir detected into dried logs without any sapwood problem on the images. For the maximum diameter measurements, they obtained a RMSE of 3.4 mm, a bias of -1.8 mm and a $R^2$ of 0.87. \emph{TEKA} run into the same order of error but works into sapwood. 

About \citet{johansson2013} algorithm, the only one in the literature dealing with knots included into sapwood, the authors presented slightly higher errors (RMSE for local diameter of 4.7 mm for pine, 5.1 for spruce) which can be due to the lower resolution of their images. Nevertheless, the errors cannot be compared directly to those of Table \ref{tabstattout} since they applied a geometric model (in the form $\phi = A + B \cdot r ^ {\frac{1}{4}}$ where $r$ is the distance to the pith and $\phi$ is the conic angle from pith to the local diameter) both to the outputs of their algorithm and to the reference measurements. With applying the same model to our data we obtain significantly lower RMSE of 3.6 mm for pine and 3.0 mm for spruce.

In the literature, some authors \citep{Breining2012, Nordmark2003, Andreu2003, Oja2000} validated their algorithm by comparison with manual measurements made on real boards or cross-sections. Instead, we decided to validate our algorithm by comparison with manual measurements on CT images \citep[like][]{Fleur01, johansson2013} and not on real wood samples. The reason of our choice is that the comparison between knot borders visible on color images (i.e., based on wood color variations) and on corresponding CT images (i.e., based on wood density variations), although very interesting, is a distinct problem, totally independent of the algorithm performance, and which should be studied separately. 
One reason for studying the correspondence between the knot borders based on color and wood density variations is that current grading rules are defined based on knot sizes measured on real boards. 

To distinguish between sound and dead knots or between sound and dead parts of knots based on CT images is not an easy task since the pattern of annual rings around knots is not always clearly visible, especially within sapwood where the moisture content is high. 
Based on our knot segmentation, another way to estimate the position of the limit would be to check the knot diameter profile radially and to put the limit where the knot has stopped its growth. However, it would not be very accurate because there is a period of time during which the knot does not grow anymore (or almost anymore) but is still sound. Some authors mention a period of decline between eight and ten years during which the knot is still alive but without growing \citep[e.g., ][]{dietrich1973, kershaw1990, Colin1992}. \cite{johansson2013} assuming that the dead knot border is located at the maximal knot diameter obtained a RMSE of 12 mm for the position of the border with a $R^2$ of 0.19 based on 65 knots of Scots pine.
In the literature, very few studies provide validation results about the detection of the dead knot border and the results are in overall not satisfactory. Moreover, \cite{johansson2013} mention that in order to accurately estimate product value, it is more valuable to improve diameter measurements than to improve the dead knot border detection.  

In addition to wet sapwood, another similar problem is the presence of wet heart which occurs very often in silver fir trees \citep{torelli2009}. The reason for the presence of wet heart in fir trees is not well known. Our sample included knots passing through wet heart (17 knots from a total of 31 fir knots) and the corresponding results were satisfactory (although not perfect) thanks to the correction of outliers. 

In this work, we made the assumption, both for manual measurements and in the algorithm, that knot cross-section was circular. Actually, it is known that knot vertical diameter is slightly bigger than knot horizontal diameter (measurements in axes perpendicular to the longitudinal axis of the knot). On Sitka spruce, \cite{achim2006} found a slight difference between horizontal and vertical diameters for whorl branches but no difference for inter-whorl branches and \citet{Merkel1967} reported a 1.057 ratio between diameter measured vertically and diameter measured horizontally for Norway spruce knots.

Five softwood species were used for the development and for the validation of the algorithm. It was thus necessary to develop a generic approach by identifying the common patterns to all these species. Another approach would have been to adapt the parameters for each species independently. Knot characteristics are very variable depending on the species : size, inclination, shape, density. We have observed such differences for the five softwood species based on a more complete sample of 1668 knots from six logs per species (results obtained from a Master thesis work not published). However, in the heartwood knots appeared denser than normal surrounding wood for all the five softwood species and in the sapwood the low density band around knots was present, with varying intensity, for all the species. 

We are not able to biologically explain the occurrence of a low density band around knots. Is this phenomenon due to a difference in moisture content or in wood density ? This question should be further studied. The visibility of this band is related to the wood density within the band, the width of the band and the resolution of the CT images. It would be interesting to decrease the resolution of images in order to find the limit above which the band remains visible.

More generally, the question of the image resolution has to be discussed. In the present paper, the images were provided by a medical CT scanner and the resolution could be considered as relatively high. Our first objective was knot measurements for scientific purposes. In an industrial context, the longitudinal resolution would probably be about 1 image every 1 cm and the transversal resolution would be in the order of 1 mm/pixel. The impact of a low longitudinal resolution will be different depending on the species due to the species-specific knot inclination and knot size. For example, some species like fir have almost horizontal and small branches whereas branches are much bigger and more inclined for pine trees. As a consequence, a knot from fir will be present on a very few number of consecutive images, often only on one image, which will limit the 3D reconstruction of knots. In the literature, \cite{johansson2013} explain that they have simulated low resolution images from images provided by a medical CT scanner in order to be comparable with the industrial context. However, they do not explain how they have decreased the resolution and what is the final transversal resolution of the images (the longitudinal resolution is 1 image every 1 cm) on which the algorithm has to work. In our eyes, a low longitudinal resolution and a lack of 3D support for knots is more limiting than a low transversal resolution. Nevertheless, the algorithm should work accurately on the biggest knots (several centimeters of diameter) which are the main concern for structural uses. We also could expect that the resolution will not be limiting anymore in the near future due to technological improvements regarding the next generations of industrial CT scanners.

The algorithm was implemented as an ImageJ plug-in programmed in Java language without any optimisation concern. It was applied here to CT images obtained with a medical scanner. Before envisaging online application in industrial conditions further work is needed:
\begin{itemize}
	\item The processing time must be reduced. It should be relatively easy since the algorithm does not require any very time consuming operation and could work in parallel on several knots;
	\item The effect of reducing image resolution on the accuracy of measurements has to be tested.
\end{itemize}

Another needed improvement would be the detection of the knot end. In the present work, all the knots were reaching the bark. It would be possible to find criteria to estimate the reliability of the knot pith and diameter detection on each transversal image. For the reliability of the pith detection, we could use the radial profile of the Hough accumulation value. A decrease in the accumulation at the bark side would indicate a potential problem in the detection of the pith or signify that the knot is no more present at this location. Future implementation of the algorithm will include this improvement.

\subsection*{Weaknesses of the proposed method}
\begin{itemize}
	\item The knot pith detection is difficult and few accurate at the stem pith side due to the low number of edge pixels (small diameter of the knot at this location) voting in the Hough accumulation method;
	\item An extrapolation is used for estimating the diameter of the last 10\% of the knot length located at the bark side;
	\item The algorithm does not detect the knot end;
	\item The algorithm depends on the knot detection step which should have previously isolated each knot within a 3D sector in view of the segmentation step;
	\item The algorithm is not yet optimized.
\end{itemize}

\subsection*{Strengths of the proposed method}
\begin{itemize}
	\item The algorithm works in the heartwood and in the sapwood;  
	\item The detection is generic and robust because it is based on the observation of a typical pattern of density variation around knots and not on a fixed threshold of density which would depend on the species and on the moisture content. Moreover the algorithm works even if the low density band is only partially visible around the knot;
	\item The algorithm is easy to optimize by parallelization of the processing of each knot.   
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
\section{Conclusion}

The knot segmentation algorithm named \textit{TEKA} was developed based on a tangential approach. 
The algorithm starts from a stack of sub-images oriented tangentially to the log annual growth rings. The stack of consecutive sub-images follows the knot radial direction.
After detecting the knot pith on each tangential slice, a polar elliptic transformation centered on the knot pith is applied for computing the profile of density from the knot pith to the normal wood surrounding the knot. The knot diameter is estimated by searching a local minimum in the profile. This minimum seems to correspond to the normal stem wood in heartwood and to a low density band that was observed around most of the knots in sapwood.

The method was applied to a set of 125 knots from five softwood species: Douglas-fir, silver fir, European larch, Scots pine and Norway spruce.
The accuracies of knot positions and diameters were assessed by comparison with manual measurements.
The errors were almost identical in heartwood and sapwood, and few differences were observed between species. Compared to previously published algorithms, \textit{TEKA} does not seem more accurate but works as well for sapwood and heartwood contrarily to other comparable algorithms. Moreover, it seems robust enough to process a large range of knot morphologies and various aspects of sapwood including more or less advanced states of drying.

A free and open source version of the \textit{TEKA} algorithm was implemented as a plug-in\footnote{\url{http://www.nancy.inra.fr/foret-bois-lerfob/Xylosciences/Tomographie-X}} for ImageJ.
A fully automated version will be soon embedded in the TKDetection software\footnote{\url{https://github.com/akrah/TKDetection}} \citep{Adrien01, Adrien02, Adrien03} which will perform the preliminary detection step and call \textit{TEKA} for each detected knot. This tool will be used in a next future for performing a sensitivity analysis on the input parameters and to validate the complete processing (i.e., knot detection and segmentation) on a larger sampling.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
\section*{Acknowledgement}
We would like to thank Ets. Siat-Braun who graciously supplied the log samples and Charline Freyburger who performed scanner measurements. The UMR 1092 LERFoB is supported by a grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (ANR-11-LABX-0002-01, Lab of Excellence ARBRE).

\section*{References}

\bibliographystyle{elsarticle-harv}
\bibliography{bibliographie}

%======================================================================

\clearpage

\begin{table*}[!h]
\footnotesize
\centering
\begin{minipage}[c]{\textwidth}
\begin{tabular}{ l | C{.7cm} 
| C{.75cm} C{.75cm} C{0.5cm} 
| C{.75cm} C{.75cm} C{0.5cm} 
| C{.75cm} C{.75cm} C{0.5cm} 
| C{.75cm} C{.75cm} C{0.5cm} 
|}

\multicolumn{2}{c}{}
& \multicolumn{3}{c}{\textbf{Maximal}}
& \multicolumn{3}{c}{\textbf{Local}}
& \multicolumn{3}{c}{\textbf{Vertical}}
& \multicolumn{3}{c}{\textbf{Horizontal}}\tabularnewline

\multicolumn{2}{c}{}
& \multicolumn{3}{c}{\textbf{diameter}}
& \multicolumn{3}{c}{\textbf{diameter}}
& \multicolumn{3}{c}{\textbf{position}}
& \multicolumn{3}{c}{\textbf{position}}\tabularnewline

\cline{2-14}
& $N$\footnote{$N$ is the number of knots. The number of observations is $N$ for maximal diameter, $N\times10$ for local diameter and $N\times11$ for vertical and horizontal positions.}
& RMSE (mm) & Mean bias (mm) & $R^2$
& RMSE (mm) & Mean bias (mm) & $R^2$
& RMSE (mm) & Mean bias (mm) & $R^2$ 
& RMSE (mm) & Mean bias (mm) & $R^2$ \tabularnewline

% Soumis (124 noeuds)
%\hline
%\multicolumn{1}{|l|}{Douglas-fir} & 4.17 & +0.61 & 0.80 & 2.57 & -0.28 & 0.99 & 0.93 & -0.01 & 0.99 \tabularnewline
%\multicolumn{1}{|l|}{Silver fir} & 2.37 & +1.20 & 0.71 & 2.22 & +0.08 & 0.92 & 1.57 & -0.09 & 0.89 \tabularnewline
%\multicolumn{1}{|l|}{European larch} & 2.15 & +0.45 & 0.92 & 6.16 & +0.81 & 0.79 & 1.43 & +0.22 & 0.96 \tabularnewline
%\multicolumn{1}{|l|}{Scots pine} & 4.19 & +1.77 & 0.86 & 5.76 & -0.47 & 0.97 & 2.20 & +0.16 & 0.89 \tabularnewline
%\multicolumn{1}{|l|}{Norway spruce} & 3.35 & -0.68 & 0.80 & 2.90 & +0.36 & 0.94 & 1.69 & -0.02 & 0.95 \tabularnewline
%\hline
%\multicolumn{1}{|l|}{All species} & 3.28 & +0.67 & 0.85 & 4.14 & +0.12 & 0.96 & 1.59 & +0.04 & 0.97 \tabularnewline
%\hline

% Révision (125 noeuds + modif de l'interpolation)
\hline
\multicolumn{1}{|l|}{Douglas-fir} & 
25 & 3.44 & +0.14 & 0.91 & 4.40 & +0.96 & 0.79 & 2.64 & -0.29 & 0.99 & 0.93 & 0.00 & 0.99 \tabularnewline%DOU
\multicolumn{1}{|l|}{Silver fir} & 
29 & 2.61 & +1.27 & 0.67 & 2.43 & +1.19 & 0.70 & 2.34 & -0.14 & 0.91 & 1.67 & -0.10 & 0.88 \tabularnewline%FIR
\multicolumn{1}{|l|}{European larch} & 
26 & 1.83 & -0.18 & 0.94 & 2.10 & +0.63 & 0.93 & 6.05 & +0.53 & 0.80 & 1.39 & +0.23 & 0.96 \tabularnewline%LAR
\multicolumn{1}{|l|}{Scots pine} & 
20 & 5.70 & +2.25 & 0.75 & 4.34 & +2.04 & 0.85 & 6.09 & -1.18 & 0.97 & 2.22 & +0.17 & 0.89 \tabularnewline%PIN
\multicolumn{1}{|l|}{Norway spruce} & 
25 & 4.91 & +0.22 & 0.64 & 3.28 & -0.56 & 0.81 & 2.77 & +0.28 & 0.94 & 1.60 & -0.02 & 0.95 \tabularnewline%SPR
\hline \multicolumn{1}{|l|}{All sapwood} & 
125 & - & - & - & 3.63 & +0.57 & 0.86 & 3.40 & +0.52 & 0.98 & 1.56 & +0.13 & 0.97 \tabularnewline%aubier
\multicolumn{1}{|l|}{All heartwood} & 
125 & - & - & - & 3.12 & +1.03 & 0.80 & 4.73 & -0.58 & 0.92 & 1.61 & -0.02 & 0.97 \tabularnewline%duramen
\hline \multicolumn{1}{|l|}{All species} & 
125 & 3.83 & +0.69 & 0.85 & 3.37 & +0.81 & 0.85 & 4.21 & -0.11 & 0.96 & 1.59 & +0.05 & 0.97 \tabularnewline%Algo
\hline
\end{tabular}
\end{minipage}
\caption{Main statistics on diameter measurements and pith positioning by species and wood compartment.}
\label{tabstattout}

\end{table*}

%======================================================================

\clearpage

\listoffigures

\clearpage

% ==== FIG 1 ====
\begin{figure*}[!h]
\centering
\includegraphics[width=14cm]{Matrice_images_billons_avec_noeuds}
\caption{Illustration of one CT slice per log in order to visualize some knots and the sapwood, as well as to obtain an information about the moisture content of logs. The three first letters of the log names identify the species, the following digit identifies the tree, the last letter indicates the stem bottom (B) or top (T).}
\label{imagesBillons}
\end{figure*}

\begin{figure*}[!h]
\centering
\includegraphics[width=\textwidth]{Images_zone_noire}
\caption{Transversal, radial and tangential sections for the five softwood species (knots DOU-1-T-58, FIR-1-T-278, LAR-4-T-96, PIN-2-T-96 and SPR-3-T-301. The knot numbers include the log name followed by a slice number).}
\label{Images_zone_noire}
\end{figure*}

% ==== FIG 3 ====
\begin{figure}[!h]
\centering
\includegraphics[width=9cm]{decoupe-tangentielle}
\caption{Illustration of the tangential reslicing method (log DOU-1-B). Parallel red lines correspond to the planes used for reslicing.}
\label{decoupeTan}
\end{figure}

% ==== FIG 4 ====
\begin{figure}[!h]
\centering
\subfloat[]
{
\includegraphics[height=5cm]{PIN-1-T-pith}
\label{polara}
}
\,
\subfloat[]
{
\includegraphics[height=5cm]{PIN-1-T-polar}
\label{polarb}
}
\,
\subfloat[]
{
\includegraphics[height=5cm]{PIN-1-T-polarellipse}
\label{polarc}
}
\caption{Example of polar transformation applied to a slice of Scots pine knot PIN-1-T-107. (a) Original tangential slice. (b) Polar circular transformation. (c) Polar elliptic transformation. The automatically detected knot pith (red cross) is the centre of the polar transformations. In the polar images, the vertical axe represents the azimuth around the pith (with a one degree step) and the horizontal axe the distance to the pith. After the elliptic transformation, the knot contour appears as a rather straight vertical line, allowing to compute the radial profile of Fig. \ref{histogramme} (sapwood) by averaging the pixel columns.
}
\label{polar}
\end{figure}

% ==== FIG 5 ====
\begin{figure}[!h]
\centering
\includegraphics[width=9cm]{histogramme}
\caption{Profile of mean Hounsfield units for two slices of Scots pine knot PIN-1-T-107. X-axis corresponds to the horizontal distance from the knot pith. Point A is the global maximum value of the profile. Point B is the first local minimum encountered after A. Point C is computed with equation \ref{eqrayon}, C abscissa giving the knot minor radius.
The sapwood slice is shown in Fig. \ref{polar} and \ref{segmentationa}, the heartwood slice is shown in Fig.  \ref{segmentationb}.
}
\label{histogramme}
\end{figure}

% ==== FIG 6 ====
\begin{figure}[!h]
\centering

\subfloat[]
{
\includegraphics[height=5cm]{PIN-1-T-seg}
\label{segmentationa}
\qquad
}
\subfloat[]
{
\includegraphics[height=5cm]{PIN-1-T-seg-duramen}
\label{segmentationb}
}
\caption{Automatic segmentation for two tangential images of Scots pine knot PIN-1-T-107. (a) Sapwood slice. (b) Heartwood slice.}
\label{segmentation}
\end{figure}

% ==== FIG 7 ====
\begin{figure}[!h]
\centering
\includegraphics[width=9cm]{lowess}
\caption{Post-processing steps applied to Scots pine knot PIN-1-B-33: (a) Raw data with two groups of outlier points; the red line is the LOWESS curve fitted on the data. (b) Outliers are removed and gaps are filled. (c) The profile is smoothed and the end of the profile is extrapolated.}
\label{lowess}
\end{figure}

% ==== FIG 8 ====
\begin{figure*}[p] %SPR-2-T-301-921 
\centering
\subfloat[89$\times$102 mm]
{
\includegraphics[width=2.3cm]{XY-150}
\includegraphics[width=2.3cm]{XY-90}
\includegraphics[width=2.3cm]{XY-60}
\includegraphics[width=2.3cm]{XY-30}
\includegraphics[width=2.3cm]{XY-01}
\label{result1a}
}
\qquad
\subfloat[89$\times$139 mm]
{
\includegraphics[width=2.3cm]{XZ-56}
\includegraphics[width=2.3cm]{XZ-64}
\includegraphics[width=2.3cm]{XZ-78}
\includegraphics[width=2.3cm]{XZ-88}
\includegraphics[width=2.3cm]{XZ-95}
\label{result1b}
}
\caption{Automatic segmentation of a 30 mm diameter knot of Norway spruce (knot SPR-2-T-301) from the log pith (left) to the bark (right). (a): Tangential views. (b): Transversal views.}
\label{result1}
\end{figure*}

% ==== FIG 9 ====
\begin{figure*}[p] % DOU-1-T-173-540
\centering

\subfloat[59$\times$71 mm]
{
\includegraphics[height=3.5cm]{DOU-1-T-tan}
\label{result2a}
}
\qquad
\subfloat[90$\times$71 mm]
{
\includegraphics[height=3.5cm]{DOU-1-T-radial}
\label{result2b}
}
\qquad
\subfloat[59$\times$90 mm]
{
\includegraphics[height=3.5cm]{DOU-1-T-transv}
\label{result2c}
}
\caption{Automatic segmentation of a 6 mm diameter knot of Douglas fir (knot DOU-1-T-173). (a) Plane AA. (b) Plane BB. (c) Plane CC. In red the extrapolated part of the knot.}
\label{result2}
\end{figure*}

% ==== FIG 10 ====
\begin{figure*}[p]
\centering
\subfloat[82$\times$163 mm]
{
\includegraphics[height=5.5cm]{PIN-1-B-33-878-tan}
\label{result3a}
}
\qquad
\subfloat[135$\times$163 mm]
{
\includegraphics[height=5.5cm]{PIN-1-B-33-878-rad}
\label{result3b}
}
\qquad
\subfloat[82$\times$135 mm]
{
\includegraphics[height=5.5cm]{PIN-1-B-33-878-trans}
\label{result3c}
}
\caption{Automatic segmentation of a 41 mm diameter knot of Scots pine (knot PIN-1-B-33). (a) Plane AA. (b) Plane BB. (c) Plane CC.  In red the extrapolated part of the knot.}
\label{result3}
\end{figure*}

% ==== FIG 11 ====
\begin{figure}[!h]
\centering
\includegraphics[width=.45\textwidth]{diams_manu-auto}
\includegraphics[width=.45\textwidth]{diamax_manu-auto}
\caption{Comparison between automatic and manual measurements of local diameters (left side) and maximal diameters (right side). Diagonal lines show the $y=x$ lines.}
\label{diameter_plot}
\end{figure}

% ==== FIG 12 ====
\begin{figure}[!h]
\centering
\includegraphics[width=9cm]{residus_diam}
\caption{Error on local diameter measurement as a function of the relative distance along the knot. The circles identify the points out of 1.5 $\times$ the interquartile range. Relative position 0 corresponds to the log pith and 100\% to the log bark.}
\label{diameter_error_boxplot}
\end{figure}

% ==== FIG 13 ====
\begin{figure}[!h]
\centering
\includegraphics[width=9cm]{moelleabsdist_40mm}
\caption{Absolute error on vertical and horizontal pith positioning as a function of the relative distance along the knot. The circles identify the points out of 1.5 $\times$ the interquartile range. Relative position 0 corresponds to the log pith and 100\% to the log bark; two outliers with vertical errors of 76 and 46 mm are not visible.}
\label{moellez}
\end{figure}

%\begin{figure*}[!h]
%\centering
%\includegraphics[width=15cm]{structure_for_profile}
%\caption{Some example where the dark profile is not needed to find the knot diameter. Arrows point the dark structures which allows to have the profile discribed in part \ref{substep2}.}
%\label{darkstructures}
%\end{figure*}


\end{document}